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Self-consistent modeling of pair cascades in the polar cap of a pulsar. 
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Here we briefly report on first results of self-consistent simulation of non- stationary electron- 
positron cascades in the polar cap of pulsar using specially developed hybrid PfC/Monte-Carlo 
numerical code. We consider the case of Ruderman-Sutherland cascade - when particles cannot be 
extracted from the surface of the neutrons star. 
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I. INTRODUCTION 

Many previously proposed models for polar cap cas- 
cades (and almost all quantitative models) assumed 
stationary particle outflow [l|, [H, d, EL HI- Predic- 
tions of such models disagree with both observational 
data (e.g. the number of electron-positron pairs in 
pulsar wind nebula is much higher than predicted, 
e.g. [6|) and results of numerical models of force- free 
pulsar magnetosphere (the current density required to 
support force-free magnetosphere differs substantially 
from what stationary model for polar cap cascade pre- 
dicts [7(). On the other hand, the stability of such sta- 
tionary models has not been quantitatively studied. 

Particle acceleration and electron-positron pair pro- 
duction in the polar cap could be essentially non- 
stationary: time intervals of effective particle accel- 
eration could alternate with intervals when the accel- 
erating electric field is screened by electron-positron 
pairs created in the cap [1, [tj [ljl H3 • To construct 
a consistent model for particle acceleration in the po- 
lar cap and high energy emission produced there we 
need to know the pattern of particle flow. In our view, 
the study of electron-positron cascades should be done 
starting ab initio. The key "ingredients" of the system 
must be included in the model: back reaction of par- 
ticles on the accelerating electric field and the delay 
between photon emission and pair injection. Possible 
complexity of the system behavior compels us to con- 
duct a numerical experiment where particle accelera- 
tion, pair production and variation in the accelerating 
electric field are modeled self-consistently. 

Here we presents first results of such self-consistent 
modeling. We consider the case of Ruderman- 
Sutherland cascade and describe the main properties 
of the discharge. 



II. NUMERICAL METHOD 

For modeling of pair cascades we developed special 
hybrid PIC/Monte Carlo code. Current version of the 
code is ID. The code flow schematic is shown in Fig. [I] 
the code works as follows. 

Plasma dynamics is calculated according to the 
standard PIC algorithm. Using the current density 
known from the previous step we solve Maxwell equa- 



tions and get the electric field at the grid points. Then 
for each particle we interpolate the electric field to the 
particle's position and get the force on the particle. By 
solving the equation of motion we advance particle's 
momentum and position. Motion of particles across 
cell boundaries is counted as their contribution to the 
electric current which is collected and stored for the 
next time step. 

Photon emission and pair production arc calculated 
as follows. We sample how many photons capable of 
producing electron-positron pairs each particle emits 
during the current time step. For each emitted pho- 
ton we sample its energy from the energy distribution 
for a given emission process. Then we sample the 
distance which the photon will travel until it gets ab- 
sorbed. Calculation of the optical depth to pair cre- 
ation is done with space steps varying according to 
the current value of the cross-section for photon ab- 
sorption. Most of the steps are much larger that the 
cell size. The energy of the photon, the position and 
the time of the absorption are stored in an array. At 
every time step we iterate over that array and pick up 
photons which must be absorbed at the time of the 
current time step. For each of the selected photons 
we inject an electron and a positron at the place of 
the absorption and delete that photon from the array. 
Being injected at the same point, the electron and the 
positron do not change charge and current densities 
at the moment of injection. 

If there are too many numerical particles of some 
kind in the computational domain, their number can 
be reduced by deleting some randomly selected par- 
ticles. Statistical weights of the selected particles are 
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FIG. 1: Flow schematic for the hybrid PIC/Monte-Carlo 
scheme. 
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summed and then the statistical weights of all remain- 
ing particles of the same kind as deleted ones are in- 
creased correspondingly in order to compensate for 
the deleted particles. 

III. PHYSICAL MODEL 

As previously there were no direct self-consistent ki- 
netic studies of time depended cascade starting from 
the first principles, we decided to address first a simple 
case in order to develop an intuition of what kind of 
plasma behavior to expect and to adjust the numerical 
technique accordingly. Ruderman and Sutherland 0] 
cascade is the simplest possible model for the pair cas- 
cade in the polar cap - there is no plasma inflow from 
the surface of the neutron star (NS) and all plasma in 
the cascade zone is produced by pair creation. 

Ruderman and Sutherland [8| estimate for the 
height of the cascade zone for young pulsars (their 
eq. (22)) gives 

h RS ~5xl0 3 pl /T P 3 / r B^ 7 cm. (1) 

For young pulsars, with periods less than ~ 0.1 sec, 
has is l ess that the width of the polar cap r pc ~ 
1.4xl0 4 /"/P cm. Therefore, one-dimensional approx- 
imation should work well for such cascades. In ID 
the continuity equation for the charge density and the 
Gauss equation for the accelerating electric field can 
be combined into the single equation for the electric 
field E (see e.g. [n|) 

9E 

_=-47r(j-jo), (2) 

where j is the actual current density and jo is the 
mean current density flowing through the system. 
This is the equation we solve for the electric field. It 
does not need explicitly set boundary conditions on E. 
Boundary conditions are implicitly set by the choice of 
jo. The initial electric field distribution is calculated 
as solution of ID Gauss equation for some given ini- 
tial charge density distribution and initial boundary 
conditions. 

For young pulsars the dominating emission pro- 
cess in terms of number of pair-production capable 
photons is the curvature radiation (e.g. @). We 
are primarily interested in the dynamics of the dis- 
charge zone (the region with the accelerating elec- 
tric field). The size of that zone should be of the 
order of few has- Synchrotron photons emitted by 
electron-positron pairs are much less energetic than 
the curvature photons and, therefore, these photons 
have much larger mean free paths. They are absorbed 
at large distances from the NS, where plasma density 
is expected to be very high and electric field is al- 
ready screened. Consecutively, pairs produced by the 
synchrotron photons do not influence the discharge 



dynamics and we ignore synchrotron emission in our 
simulations. 

So, our model for the Ruderman-Sutherland cas- 
cade includes ID electrodynamics, curvature radiation 
as the photon emission process, photon absorption 
and pair creation in strong magnetic field. Particle 
equation of motion includes radiation reaction due to 
curvature radiation. 



IV. RESULTS OF NUMERICAL 
SIMULATIONS 

Numerical simulation have shown that in the 
Ruderman-Sutherland model pair creation is quasi- 
periodic and self-sustained. We performed simulations 
for different initial particle distributions and different 
initial electric fields, different strengths of the mag- 
netic field, different radii of curvature of magnetic field 
lines and for different pulsar periods. After the initial 
burst of pair creation the cascade zone always settled 
down to a quasi-periodic behavior qualitatively simi- 
lar for all physical parameters admitting pair creation. 
The system seems to forget initial conditions and af- 
ter a short relaxation (couple of flyby times) for given 
magnetic field, pulsar period and the mean current 
jo its behavior is the same independent on the initial 
configuration. 

We describe here main properties of the cascade us- 
ing as an example the case with the mean current den- 
sity equal to the Goldreich- Julian [l2j current density, 
jo = jcj ■ For the mean current density different from 
j G] cascade behavior is qualitatively similar. In this 
example we consider a pulsar with period P — 0.2 sec, 
magnetic field B = 10 12 G, radius of curvature of the 
magnetic field lines pcur = 10 6 cm (such value for pcur 
was used in [8j). We performed simulations for pure 
dipole magnetic field with pcur ~ 10 8 cm too. Quali- 
tatively results do not depend on the radius of curva- 
ture, but for smaller pcur calculations with the same 
numerical resolution can be done faster, as the size 
of the gap with accelerating electric field is smaller. 
Angular velocity of NS rotation is anti-parallel to the 
magnetic moment of the star, so the Goldreich-Julian 
change density is positive. 

We describe below a whole cycle of pair formation 
and illustrate the cascade development by a series of 
snapshots at several time moments during one cy- 
cle with plots showing different physical quantities in 
Figs. [31 |H In these figures we present plots for a 
typical cycle of pair formation taken from a long sim- 
ulation where several such cycles were observed. The 
time on this figures is normalized to the flyby time of 
the computational domain. Time is counted from the 
start of a particular simulation, so its absolute value 
has no physical meaning - only time intervals between 
the shots have physical meaning. 

Changes in the charge density distribution gives the 



eConf C091122 



2009 Fermi Symposium, Washington, D.C., Nov. 2-5 



3 
2 

P 1 

-1 



-jt 



-30 - 



Jt= 10.100 L , , , 


"ill 




I 





.1 .2 .1 .2 .1 .2 .1 .2 .1 .2 




FIG. 2: Snapshots of a full pair-formation cycle. Charge density is normalized to the Goldreich-Julian charge density. 
Distance is normalized to the polar cap radius r pc , NS surface is on the left, at % = 0. Pulsar period P — 0.2 sec, 
magnetic field B = 10 12 G, radius of curvature of the magnetic field lines pcur = 10 6 cm. Time, shown in a small box 
in the upper left corner of each plot, is normalized to the flyby time of the computation domain. The presented cycle is 
taken from the middle of a long simulation. A typical discharge is shown, top The whole cycle of cascade development 
is shown at equally separated moments of time; bottom The formation of pair plasma blob, marked by the gray area 
on the top panel, is shown here with smaller time intervals between snapshots. 



best overview of what is going on in the discharge 
zone, because charge density indirectly provides in- 
formation about both the particle number density and 
the electric field. In Fig. [2] we plot the change density 
at equally spaced time interval during the discharge 
cycle. In the upper panel of that figure we present 
an overview of the entire cycle, in the lower panel we 
plot snapshots of the change density distribution at 
smaller time intervals for the most interesting part of 



the discharge - formation of a new plasma blob. In 
Figs. [3l 3] we show more detailed information about 
physical conditions in the discharge zone: the number 
densities of electrons and positron, the accelerating 
electric field, phase portraits of electrons, positrons 
and pair producing photons. On the phase portraits 
particles with positive values of the 4- momentum p are 
those which move from the NS, particles with negative 
p move toward the NS. 
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FIG. 3: Snapshots of cascade development at time steps marked on the bottom panel of Fig. [2] by the yellow time boxes. 
At each moment of time several characteristic of the cascade zone are shown as functions of the distance from the NS: 
I s row: charge density of electrons (blue) and positrons (red), p e ±, normalized to the Goldreich- Julian charge density; 
2 nd row: total charge density p normalized to the Goldreich- Julian charge density; 3 rd row: accelerating electric field 
Eu normalized to the vacuum electric field; 4 th row: Phase space of positrons (normalized to m e c momentum p e + vs 
coordinate x); 5 th row: Phase space of electrons (normalized to m e c momentum p e - vs coordinate x); 6 th row: Phase 
space of pair-producing gamma-rays (normalized to m e c momentum p 1 vs coordinate x). 



A typical cycle of the discharge starts with a vac- 
uum gap forming above the surface of the NS (snap- 
shots with t — 8.9 — 9.5). The electric field there 
is very strong and charged particles entering the gap 
are accelerated up to very high energies. Plasma cre- 
ation starts close to the NS and is ignited by the 
gamma-rays emitted by electrons flowing toward the 
NS. These "primary" electrons have been created in 
the previous bursts of pair formation. They leaks from 
the tail of the plasma blob created in the previous cy- 



cle and enters the gap from above. The newly cre- 
ated electrons and positrons are accelerated by the 
strong electric field of the gap and start producing 
high energy photons, many of which decay into pairs 
too (snapshot with t = 9.667 in Figs. [2]). In Fig.rjH 
on the plots showing the phase portrait for electrons 
"primary" electrons accelerated in the gap can be seen 
as the particle population having thin line-like form. 
Particle populations scattered over the large area be- 
ginning at the left end of the phase space plots repre- 
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FIG. 4: Snapshots of cascade development at time steps marked on the bottom panel of Fig. [2] by the yellow time boxes 
(continued from Fig. [3}. Notations are the same as in Fig. [3] Note the change of y-axis normalization for charge densities 
(the first two rows) compared to Fig. El 



sent secondary particles. 

While number density of the pair plasma remains 
less then the Goldreich- Julian density, the electric 
field remains strong and electrons and positrons are 
accelerated to energies high enough to emit photons 
capable of producing pairs (see snapshots with t = 
9.667 — 9.8). When the number density become larger 
than the Goldreich- Julian density, acceleration of par- 
ticles ceases and particles created after that moment 
do not emit pair-producing photons. Screening of the 
electric field begins near the NS surface, where first 
pairs are formed, so the region of a very low plasma 
density (the gap) with strong electric field detaches 
from the NS surface and propagates into the magne- 



tosphere (snapshots with t > 9.933). From below this 
gap is limited by the blob of freshly formed plasma, 
from above - by the plasma created in the previous 
burst of pair formation. Particles from the previ- 
ous bursts of pair formation which have momentum 
directed toward the NS enters the gap from above. 
Electrons are accelerated toward the NS, positrons are 
turned back by the strong electric field of the gap and 
move into the magnetosphere. Because of this rever- 
sal the gap upper boundary moves noticeably slower 
than the speed of light. The front of the blob with 
the freshly created plasma consists of ultra-relativistic 
positrons accelerated in then strong electric field, so it 
moves relativistically and, therefore, the gap shrinks 
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when the blob moves into the magnetosphere. Pairs 
are continuously injected into the blob, because it 
practically co-moves with the pair-producing photons. 
Eventually the front of the new blob catches the tail of 
the previously created blob. Therefore, the magneto- 
sphere will be filled with plasma and there will be no 
gaps in plasma spacial distribution far from the polar 
cap. When the blob has traveled some distance a new 
gap starts forming, ignited by the electrons leaving 
the current blob. 

When the electric field is still strong large amplitude 
plasma oscillations are excited in the newly formed 
plasma. These oscillations persists after the global 
accelerating field has been screened. Fluctuating elec- 
tric field of the oscillations reverses some electrons and 
positrons toward the NS. These particles forms the 
above mentioned tail of the plasma blob. The elec- 
trons from that tail will be the "primary" particles in 
the next burst of pair formation. 

The height of the gap is ~ 2 times larger than esti- 
mate given by eq. |T]). The maximum particle energy 
is ~ 4 higher than given in Q and particle energy 
distributions is broad. 



V. CONCLUSIONS AND OPEN ISSUES 

Main findings from the numerical simulations can 
be summarized as follows. The Ruderman-Sutherland 
cascade can operate for any positive current density. 
The cascade is self-sustained and discharges occur 
quasi-periodically, the whole systems shows a sort of 
limit cycle behavior. As cascade properties do not de- 
pend on the initial conditions in our simulations, they 
are uniquely set by the mean current density jo, the 
potential drop across the polar cap and the curvature 
of the magnetic field lines. Distribution of the pair 
plasma in the magnetosphere is continuous, with no 
gaps, so the parallel component of the electric field far 
from the polar cap should be screened. Thus, such cas- 
cade model agrees well with force-free magnetosphere 



models [e.g. 0, [H, 03 • Strong plasma oscillations are 
excited in the plasma blob during its formation. This 
might have implication for generation of radioemis- 
sion. 

While the general structure of the flow is evident 
from the performed simulations, some questions re- 
main unanswered. The most important one is about 
the period of these discharge. It depends on the rate of 
plasma leakage from the blob, i.e. how many particles 
are reversed toward the NS. The more particles leaks 
out, the later the next gap forms. Due to continu- 
ous pair injection plasma density in the blob increases 
enormously, and at some time the numerical scheme 
stops resolving the Debye length of the plasma, and 
the damping of the plasma oscillations cannot be cal- 
culated reliably. At that time results start depending 
on the numerical resolution. Because of this the blob 
cannot be followed for time interval long enough (and 
traveled distances large enough) to get the repetition 
rate of the cascade. In the presented simulations the 
size of the simulation zone is set such that the blob 
leaves calculation domain before the numerical scheme 
fails to model it correctly. So, while particle energy 
distribution can be inferred from current simulations, 
pair multiplicity and particle fluxes are known up to 
some numerical factor which depends on the repeti- 
tion rate of pair creation bursts. 

More detailed description of the results as well as 
discussion of their implication for pulsar physics will 
be given in a subsequent publication. 
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